% matlab example of VAR
% > Econometric Toolbox > Model Selection > Specification Testing

% 1. load data
% clc; clear all;
% [val txt] = xlsread('../data/data_ls_macro.xlsx',2); save '../data/data2';
% fprintf('saved data');
%%
%
clc; clear all;load '../data/data2';
GDP = val(:,3); %real GDP
M2 = val(:,4);
Credit = val(:,5);
Deposit = val(:,6);
CPI = val(:,8);
lsLNH = val(:,2);
dates = datenum(txt(2:end,1));

%%

Y = [Credit./CPI,Deposit./CPI, CPI, lsLNH];
% X = dates;

% Transforming Data for Stationarity
% plot the data for trend
figure
subplot(2,2,1)
plot(dates,Y(:,1),'r');
title('Credit')
datetick('x')
grid on
subplot(2,2,2);
plot(dates,Y(:,2),'b');
title('Deposit')
datetick('x')
grid on
subplot(2,2,3);
plot(dates, Y(:,3), 'k')
title('CPI')
datetick('x')
grid on
subplot(2,2,4);
plot(dates, Y(:,4), 'magenta')
title('interbank interest rate 1M')
datetick('x')
grid on
hold off


% transfromed and plot again
Y = [diff(log(Y(:,1:2))), diff(Y(:,3)) diff(log(Y(1:end,4)))];
X = dates(2:end);

figure
subplot(2,2,1)
plot(X,Y(:,1),'r');
title('Credit')
datetick('x')
grid on
subplot(2,2,2);
plot(X,Y(:,2),'b');
title('Deposit')
datetick('x')
grid on
subplot(2,2,3);
plot(X, Y(:,3),'k'),
title('CPI')
datetick('x')
grid on
subplot(2,2,4);
plot(X, Y(:,4),'magenta'),
title('interbank 1M')
datetick('x')
grid on


% plot in one figure only

Y(:,1:2) = 100*Y(:,1:2);
figure
plot(X,Y(:,1),'r');
hold on
plot(X,Y(:,2),'b');
datetick('x')
grid on
plot(X,Y(:,3),'k');
plot(X,Y(:,4)*10,'magenta');
legend('Credit','Deposit','CPI','interbank 1M');
hold off

%% check stationary
% temp = diff(log(lsLNH));
temp = diff(log(Deposit./CPI));
% temp = diff(CPI);
fdates = size(dates,1) - size(temp,1);
x = dates(fdates+1:end);


close all;plot(x,temp);datetick('x',27);axis tight;

[x1, pvalue] = adftest(temp); [x1, pvalue]
% % -> result :[ lsLNH CPI, log(Credit/CPI), log(Deposit/CPI) ] -> I(1)
% %           : GDP (adjusted and stationary) -> I(2)
% 
%% === SELECTING MODEL
dCredit = 100*diff(log(Credit./CPI));
dDeposit = 100*diff(log(Deposit./CPI));
dCPI = diff(CPI);
% dlsLNH = diff(lsLNH);
dlsLNH = 10*diff(log(lsLNH));
% Y = [dCredit dDeposit dCPI dlsLNH];
Y = [dCredit dDeposit dlsLNH];
numPredictors = size(Y,2);


% %%create 4 model
% dt = logical(eye(4));
% VAR2diag = vgxset('ARsolve',repmat({dt},2,1),...
%     'asolve',true(4,1),'Series',{'Credit','Deposit','CPI','1-mo interbank'},'nAR',2);
% VAR2full = vgxset(VAR2diag,'ARsolve',[]);
% VAR4diag = vgxset(VAR2diag,'nAR',4,'ARsolve',repmat({dt},4,1));
% VAR4full = vgxset(VAR2full,'nAR',4);
% VAR3full = vgxset(VAR2full,'nAR',3);
% 
% % === Choosing Presample, Estimation, and Forecast Periods. 
% 
% Ypre = Y(1:4,:);
% T = floor(.9*size(Y,1));
% Yest = Y(5:T,:);
% YF = Y((T+1):end,:);
% TF = size(YF,1);
% 
% % fitting
% [EstSpec1,EstStdErrors1,LLF1,W1] = ...
%     vgxvarx(VAR2diag,Yest,[],Ypre,'CovarType','Diagonal');
% [EstSpec2,EstStdErrors2,LLF2,W2] = ...
%     vgxvarx(VAR2full,Yest,[],Ypre);
% [EstSpec3,EstStdErrors3,LLF3,W3] = ...
%     vgxvarx(VAR4diag,Yest,[],Ypre,'CovarType','Diagonal');
% [EstSpec4,EstStdErrors4,LLF4,W4] = ...
%     vgxvarx(VAR4full,Yest,[],Ypre);
% [EstSpec5,EstStdErrors5,LLF5,W5] = ...
%     vgxvarx(VAR3full,Yest,[],Ypre);

% %% Checking Model Adequacy
% 
% [isStable1,isInvertible1] = vgxqual(EstSpec1);
% [isStable2,isInvertible2] = vgxqual(EstSpec2);
% [isStable3,isInvertible3] = vgxqual(EstSpec3);
% [isStable4,isInvertible4] = vgxqual(EstSpec4);
% [isStable5,isInvertible5] = vgxqual(EstSpec5);
% [isStable1,isStable2,isStable3,isStable4, isStable5]
% 
% % Likelihood Ratio Tests
% % count the active number parameter
% [n1,n1p] = vgxcount(EstSpec1);
% [n2,n2p] = vgxcount(EstSpec2);
% [n3,n3p] = vgxcount(EstSpec3);
% [n4,n4p] = vgxcount(EstSpec4);
% [n5,n5p] = vgxcount(EstSpec5);
% % test
% reject1 = lratiotest(LLF2,LLF1,n2p - n1p)
% reject3 = lratiotest(LLF4,LLF3,n4p - n3p)
% reject4 = lratiotest(LLF4,LLF2,n4p - n2p)
% reject5 = lratiotest(LLF5,LLF2,n5p - n2p)
% 
% 
% % AIC 
% AIC = aicbic([LLF1 LLF2 LLF3 LLF4 LLF5],[n1p n2p n3p n4p n5p])
% 
% %%
% % Compare forecast with forecast period data
% Yest = Y(5:T,:);
% [FY1,FYCov1] = vgxpred(EstSpec1,TF,[],Yest);
% [FY2,FYCov2] = vgxpred(EstSpec2,TF,[],Yest);
% [FY3,FYCov3] = vgxpred(EstSpec3,TF,[],Yest);
% [FY4,FYCov4] = vgxpred(EstSpec4,TF,[],Yest);
% [FY5,FYCov5] = vgxpred(EstSpec5,TF,[],Yest);
% 
% %% calculate error: rolling forecasting
% numPredictors = size(Y,2);
% FY1 = zeros(TF,numPredictors);
% FY2 = zeros(TF,numPredictors);
% FY3 = zeros(TF,numPredictors);
% FY4 = zeros(TF,numPredictors);
% FY5 = zeros(TF,numPredictors);
% Yest = Y(5:T,:);
% for i = 1:TF
%     FY1(i,:) = vgxpred(EstSpec1,1,[],Yest);
%     FY2(i,:) = vgxpred(EstSpec2,1,[],Yest);
%     FY3(i,:) = vgxpred(EstSpec3,1,[],Yest);
%     FY4(i,:) = vgxpred(EstSpec4,1,[],Yest);
%     FY5(i,:) = vgxpred(EstSpec5,1,[],Yest);
%     
%     if(i~=TF)
%         Yest = Y(5:T+i,:);
%         %estimate again
%         [EstSpec1,EstStdErrors1,LLF1,W1] = ...
%             vgxvarx(VAR2diag,Yest,[],Ypre,'CovarType','Diagonal');
%         [EstSpec2,EstStdErrors2,LLF2,W2] = ...
%             vgxvarx(VAR2full,Yest,[],Ypre);
%         [EstSpec3,EstStdErrors3,LLF3,W3] = ...
%             vgxvarx(VAR4diag,Yest,[],Ypre,'CovarType','Diagonal');
%         [EstSpec4,EstStdErrors4,LLF4,W4] = ...
%             vgxvarx(VAR4full,Yest,[],Ypre);
%         [EstSpec5,EstStdErrors5,LLF5,W5] = ...
%             vgxvarx(VAR3full,Yest,[],Ypre);
%     end
% end

% %% calculate error
% 
% figure
% % vgxplot(EstSpec2,Yest,FY2,FYCov2)
% % vgxplot(EstSpec4,Yest,FY4,FYCov4)
% % Sum Square Error
% error1 = YF - FY1;
% error2 = YF - FY2;
% error3 = YF - FY3;
% error4 = YF - FY4;
% error5 = YF - FY5;
% 
% SSerror1 = error1(:)' * error1(:);
% SSerror2 = error2(:)' * error2(:);
% SSerror3 = error3(:)' * error3(:);
% SSerror4 = error4(:)' * error4(:);
% SSerror5 = error5(:)' * error5(:);
% figure
% % bar([SSerror1 SSerror2 SSerror3 SSerror4 SSerror5],.5)
% bar([SSerror2 SSerror4 SSerror5],.5)
% ylabel('Sum of squared errors')
% set(gca,'XTickLabel',...
%     {'AR2 full' 'AR4 full','AR3 full'})
% title('Sum of Squared Forecast Errors')
% 
% vgxdisp(EstSpec2)

%% SELECTING MODEL: choosing lag between 1-12:
Ybackup = Y;
% Y = Y(:,[1,2,4]);

T = floor(.9*size(Y,1));

Ypre1 = Y(1:12,:);
Yest1 = Y(13:T,:);
result = zeros(12,4); %likelihood and number active parameters
YF = Y((T+1):end,:);
TF = size(YF,1);
numPredictors = size(Y,2);


for iLag = 1: 12
    VARtemp = vgxset('asolve',true(numPredictors,1),'nAR',iLag);
%     VARtemp = vgxset(VARtemp,'Series',{'Credit','Deposit','1-mo interbank'});
    [EstSpecTemp,EstStdErrorsTemp,LLFTemp,WTemp] = vgxvarx(VARtemp,Yest1,[],Ypre1);
    [isStable,isInvertible] = vgxqual(EstSpecTemp);
    [nTemp,nTempp] = vgxcount(EstSpecTemp);
    result(iLag,1:3) = [LLFTemp, nTempp,isStable];
    
    % calculate rolling forecast performance
    

    FY1 = zeros(TF,numPredictors);
    
    numPredictors = size(Y,2);
    EstSpecTemp1 = EstSpecTemp;
    Yest2 = Y(13:T,:);
    for i = 1:TF
        FY1(i,:) = vgxpred(EstSpecTemp1,1,[],Yest2);        
        if(i~=TF)
            Yest2 = Y(13:T+i,:);
%             estimate again
            [EstSpecTemp1,EstStdErrorsTemp1,LLFTemp1,WTemp1] = ...
                vgxvarx(VARtemp,Yest2,[],Ypre1);
        end
    end
    error1 = YF - FY1;
    result(iLag,4) = error1(:)' * error1(:);
end
[AIC,BIC] = aicbic(result(:,1),result(:,2),size(Yest1,1));
HQC = 2.*log(log(size(Yest1,1))).*result(:,2) - 2.*result(:,1);
xxx = [AIC BIC HQC result(:,4)];


close all;figure
plot(1:12,xxx(:,1:3)); hold on;
bar(1:12, xxx(:,4))
xlabel('Lag')
title('Information Criterion for VAR''s lag selection')
legend('AIC','BIC','HQC','error');

%% Forecasting: create best model
bestLag = 1;
finalModelSpec = vgxset('asolve',true(numPredictors,1),'nAR',bestLag);
finalModelSpec = vgxset(finalModelSpec,'Series',{'Credit','Deposit','1-mo interbank'});
[EstSpecFinal,EstStdErrorsFinal,LLFFinal,WFinal] = vgxvarx(finalModelSpec,Y(bestLag+1:end,:),[],Y(1:bestLag,:));
vgxdisp(EstSpecFinal,EstStdErrorsFinal);

%%
bestLag = 1;
Nforecast = 10;
finalModelSpec = vgxset('asolve',true(numPredictors,1),'nAR',bestLag);
finalModelSpec = vgxset(finalModelSpec,'Series',{'Credit','Deposit','1-mo interbank'});

% vgxdisp(EstSpecFinal,EstStdErrorsFinal);

ypred = zeros(Nforecast, numPredictors);
for i = 1:Nforecast
    [EstSpecFinal,EstStdErrorsFinal,LLFFinal,WFinal] = vgxvarx(finalModelSpec,Y(bestLag+1:end+i-Nforecast-1,:),[],Y(1:bestLag,:));
    [ypred(i, :), ycov] = vgxpred(EstSpecFinal,1,[],Y(1:end+i-Nforecast-1,:));    
end



% [ypred, ycov] = vgxpred(EstSpecFinal,Nforecast,[],Y(1:end-Nforecast-1,:));
y_raw = [Credit./CPI, Deposit./CPI, lsLNH];
yfirst = [log(y_raw(end-Nforecast-1,1:2)) ,log(y_raw(end-Nforecast-1,3))];
ypred1 = [yfirst;ypred(:,1:2)/100 ypred(:,3)/10];
ypred2 = cumsum(ypred1);
ypred2(:,[1,2,3]) = exp(ypred2(:,[1,2,3]));
% lasttime = dates(end);
% % timess = lasttime:31:lasttime+310;

figure;
for i = 1:3
    subplot(2,2,i);
    plot(69:79,y_raw(69:79,i),'b.'); hold on
    plot(69:79,ypred2(:,i),'r-');hold off
    title(finalModelSpec.Series{i});    
    legend('Data','Forecast','Location','NorthWest');
% %     draw vertical line
%     hx = graph2d.constantline(lasttime, 'LineStyle',':', 'Color',[.7 .7 .7]);
%     changedependvar(hx,'x');
%     datetick('x','mm-yy');
end


%% forecast with no error bars

[ypred, ycov] = vgxpred(EstSpecFinal,10,[],Y(end-bestLag-1:end,:));
y_raw = [Credit./CPI, Deposit./CPI, CPI, lsLNH];
yfirst = [log(y_raw(end,1:2)) ,y_raw(end,3), log(y_raw(end,4))];
ypred1 = [yfirst;ypred(:,1:2)/100 ypred(:,3) ypred(:,4)/10];
ypred2 = cumsum(ypred1);
ypred2(:,[1,2,4]) = exp(ypred2(:,[1,2,4]));
lasttime = dates(end);
timess = lasttime:31:lasttime+310;

figure;
for i = 1:4
    subplot(2,2,i);
    plot(dates,y_raw(:,i)); hold on
    plot(timess,ypred2(:,i),'r-');hold off
    title(finalModelSpec.Series{i});    
    legend('Data','Forecast','Location','NorthWest');
% %     draw vertical line
    hx = graph2d.constantline(lasttime, 'LineStyle',':', 'Color',[.7 .7 .7]);
    changedependvar(hx,'x');
    datetick('x','mm-yy');
end


% forecasting with simulation
rng(1);
Nsimumation = 10000;
ysim0 = vgxsim(EstSpecFinal,10,[],YF,[],Nsimumation);

y_raw = [Credit./CPI, Deposit./CPI, CPI, lsLNH];
yfirst = [log(y_raw(end,1:2)) ,y_raw(end,3),log(y_raw(end,4))];
ysim0(:,1:2,:) = ysim0(:,1:2,:)./100;
ysim0(:,4,:) = ysim0(:,4,:)./10;

ysim = [repmat(yfirst,[1,1,Nsimumation]);ysim0];
ysim2 = cumsum(ysim(:,1:end,:));
ysim2(:,[1, 2, 4],:) = exp(ysim2(:,[1,2,4],:));

ymean = mean(ysim2,3);
ystd = std(ysim2,0,3);


figure;
for i = 1:4
    subplot(2,2,i);
    plot(dates,y_raw(:,i)); hold on
    plot(timess,ymean(:,i),'r-');
    plot(timess,ymean(:,i)-ystd(:,i),'g--');
    plot(timess,ymean(:,i)+ystd(:,i),'g--');
    hold off
    title(finalModelSpec.Series{i});    
    legend('Data','Forecast','Location','NorthWest');
% %     draw vertical line
    hx = graph2d.constantline(lasttime, 'LineStyle',':', 'Color',[.7 .7 .7]);
    changedependvar(hx,'x');
    datetick('x','mm-yy');
end

%% Calculating Impulse Responses
rng(1);
fdates = arrayfun(@(x)addtodate(dates(end),x,'month'),1:20)';

FT = numel(fdates);
W0 = zeros(FT,numPredictors);
W1 = zeros(FT,numPredictors,numPredictors);
resDif = zeros(FT,numPredictors.*numPredictors);

for i = 1 : numPredictors    
    W1(1,i,i) = sqrt(EstSpecFinal.Q(i,i));
    Yimpulse = vgxproc(EstSpecFinal,W1(:,:,i),[],Y);
    Ynoimp = vgxproc(EstSpecFinal,W0,[],Y);
    
    Yimp1 = [exp(cumsum(Yimpulse(:,1:2)/100)), cumsum(Yimpulse(:,3)) ...
                exp(cumsum(Yimpulse(:,4)/10))];
    Ynoimp1 = [exp(cumsum(Ynoimp(:,1:2)/100)), cumsum(Ynoimp(:,3)) ...
                exp(cumsum(Ynoimp(:,4)/10))];
    RelDiff = (Yimp1 - Ynoimp1)./Yimp1;
    
    
    resDif(:,(i-1)*numPredictors+1:i*numPredictors) = RelDiff;
    
    figure(i);        
    subplot(2,2,1);plot(1:20,100*RelDiff(:,1));dateaxis('x',12);
    title(sprintf('Shock on %s to %s',EstSpecFinal.Series{i},EstSpecFinal.Series{1}));
    ylabel('%Change')
    
    subplot(2,2,2);plot(1:20,100*RelDiff(:,2));dateaxis('x',12);
    title(sprintf('Shock on %s to %s',EstSpecFinal.Series{i},EstSpecFinal.Series{2}))
    ylabel('%Change')
    subplot(2,2,3);plot(1:20,100*RelDiff(:,3));dateaxis('x',12);
    title(sprintf('Shock on %s to %s',EstSpecFinal.Series{i},EstSpecFinal.Series{3}))
    ylabel('%Change')
    subplot(2,2,4);plot(1:20,100*RelDiff(:,4));dateaxis('x',12);
    title(sprintf('Shock on %s to %s',EstSpecFinal.Series{i},EstSpecFinal.Series{4}))
    ylabel('%Change')
end
%%
figure(numPredictors+1);
for i = 1: numPredictors
    subplot(2,2,i);plot(1:20,resDif(:,numPredictors*i)*100);
    title(sprintf('Shock on %s to interbank rate',EstSpecFinal.Series{i}))
    ylabel('%Change')
end